import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500))
import numpy as np
import scipy.signal

4. Efectos del muestreo y fenómeno de aliasing

4.1. El espectro de una señal discreta

4.1.1. ¿Qué le ocurre al espectro de una señal continua cuando la muestreamos?

Muestrear es equivalente a multiplicar una señal continua por un tren de impulsos, también conocido como “peineta de Dirac”. La definición matemática del tren de impulsos es

\[ \upuparrows(t) = \sum_{m=-\infty}^\infty \delta[t - m / F_s], \]

donde \(F_s\) es la frecuencia de muestreo, es decir el recíproco de \(T_s\).

La siguiente figura representa graficamente el proceso de muestreo:

../../_images/dirac-comb.png

que matemáticamente sería

\[ s(t) ~ \cdot \upuparrows(t) = \{s_m\} \]

donde \(\{s_m\}\) es la secuencia muestreada a partir de \(s(t)\).

Nos interesa entender que ocurre con \(S(f)\), el espectro de \(s(t)\), cuando se realiza el muestreo. Para verificar esto matematicamente necesitamos conocer la transformada de Fourier del tren de impulsos, que es

\[ \mathbb{FT}[\upuparrows(t)] = F_s \sum_{m=-\infty}^\infty \delta(f - m F_s) \]

es decir otro tren de impulsos pero en frecuencia en lugar de tiempo. Pueden revisar la demostración de la transformada anterior aquí

Con esto tenemos todos los ingredientes para calcular la transformada de Fourier de la señal discreta

\[\begin{split} \begin{align} \mathbb{FT}[s(t) \cdot \upuparrows(t)] &= \mathbb{FT}[s(t)] * \mathbb{FT}[\upuparrows(t)] \nonumber \\ &= S(f) * F_s \sum_{m=-\infty}^\infty \delta(f - m F_s) \nonumber \\ &= F_s \sum_{m = -\infty}^{\infty} S(f - m F_s) \nonumber \end{align} \end{split}\]

Nota

El espectro de la señal discreta es equivalente a una repetición infinita del espectro de la señal continua. Las repeticiones del espectro están separadas en \(F_s\) [Hz]

Para obtener el espectro de la señal discreta hemos usado una importante propiedad de la transformada de Fourier: La transformada de Fourier de una multiplicación de señales corresponde a la convolución entre las transformadas de Fourier de cada uno de los términos.

A continuación veremos en mayor detalle que es la operación de convolución y que significa convolucionar con un tren de impulsos

4.1.2. Convolucón con un tren de impulsos

La convolución entre dos señales discretas (de una dimensión) se define matematicamente como

\[ (f * g)[n] = \sum_{m=-\infty}^\infty f[m] g[n-m] = \sum_{m=-\infty}^\infty f[n-m] g[m] \]

La variable \(m\) se suele llamar (lag). La convolución consiste en sumar la multiplicación de \(f\) con versiones retrasadas de \(g\) (o viceverza). El resultado de la convolución es una nueva función

En el caso de que \(f\) sea un tren de impulsos ocurrirá una “repetición” de \(g\) (o viceverza). La siguiente animación ejemplifica lo anterior con la convolución entre una señal gaussiana y un tren de impulsos con \(F_s=2\) [Hz]

def señal(t, mu=0, sigma=0.2):
    return np.exp(-0.5*((t-mu)/sigma)**2)

def peineta(t, Fs=2000):
    s = np.zeros_like(t)
    s[::Fs] = 1 
    return s
    
x = np.arange(-2, 2+1e-3, step=1e-3)
conv_s = np.convolve(señal(x), peineta(x), mode='same')
%%output holomap='gif' fps=5

hmap1 = hv.HoloMap(kdims=['frame'])
hmap2 = hv.HoloMap(kdims=['frame'])
tren_plot = hv.Curve((x, peineta(x)), label='Tren de impulsos')
conv_plot = hv.Curve((x, conv_s), 'Tiempo [s]', 'Convolución')
for frame in range(0, 4001, 250):    
    hmap1[frame] = hv.Curve((x, señal(x, mu=x[frame])), 'Tiempo [s]', 'Señal', label='Gaussiana')
    hmap2[frame] = hv.Scatter((x[frame], conv_s[frame])).opts(color='k', size=10)
    
main_plot = ((hmap1 * tren_plot) + (hmap2 * conv_plot)).cols(1)
main_plot.opts(hv.opts.Curve(height=200, xlim=(-2.1,2.1)))

Nota

Cada repetición de la gaussiana calza con uno de los impulsos del tren de impulsos

4.1.3. Ejemplo: Espectro de una señal Gaussiana discreta, parte 1

La señal Gaussiana con media cero es

\[ s(t) = \exp \left ( -\frac{t^2}{2\sigma^2} \right) \]

y su transformada de Fourier es

\[ S(f) = \mathbb{FT}[s(t)] = \sqrt{2\pi}\sigma \cdot \exp \left ( -2\sigma^2 \pi^2 f^2 \right). \]

Puedes ver la demostración de esta transformada aquí

Si prestamos atención a las fórmulas podemos notar que el espectro de una gaussiana es otra gaussiana pero con ancho inversamente proporcional a la original. Es decir que mientras más “ancha” sea la gaussiana en el tiempo (\(\sigma\) pequeño) más “angosto” será su espectro en frecuencia, como muestra la siguiente gráfica

x = np.arange(-5, 5, step=0.001)
gauss_time = lambda time, sigma: np.exp(-0.5*(time/sigma)**2)
gauss_freq = lambda freq, sigma: np.exp(-2*(np.pi*freq*sigma)**2)#*(np.sqrt(2*np.pi)*sigma)
gt = hv.HoloMap(kdims=['sigma'])
gf = hv.HoloMap(kdims=['sigma'])

for sigma in [0.25, 1/np.sqrt(2*np.pi), 1]:
    gt[sigma] = hv.Curve((x, gauss_time(x, sigma)), 'Tiempo [s]', 's(t)')
    gf[sigma] = hv.Curve((x, gauss_freq(x, sigma)), 'Frecuencia [Hz]', 'S(f)')
    
(gt + gf).opts(hv.opts.Curve(width=300))

Como vimos anteriormente el espectro de una señal discreta es el espectro de la señal continua convolucionado con un tren de impulsos, donde los impulsos están separados entre sí por \(F_s\) [Hz].

La siguiente figura interactiva muestra la convolución del espectro teórico de la señal gaussiana con un tren de impulsos de \(F_s=2\) [Hz] para distintos valores de \(\sigma\).

def periodic_gauss_freq(f, sigma, Fs=2):
    # Equivalente a convolucionar la gaussiana con un tren de impulsos
    S = np.zeros_like(f)
    for m in range(-20, 20):
        S += gauss_freq(f - Fs*m, sigma)
    return S

x = np.arange(-5, 5, step=0.001)
espectro_continuo = hv.HoloMap(kdims=['sigma'], sort=False)
espectro_discreto = hv.HoloMap(kdims=['sigma'], sort=False)

for sigma in [1, 1/np.sqrt(2*np.pi), 0.25]:
    espectro_continuo[sigma] = hv.Curve((x, gauss_freq(x, sigma)), 
                                        'Frecuencia [Hz]', 'Espectro', label='Continuo').opts(line_width=4)
    espectro_discreto[sigma] = hv.Curve((x, periodic_gauss_freq(x, sigma)), 
                                        'Frecuencia [Hz]', label='Discreto')
    
(espectro_continuo * espectro_discreto * hv.Box(0, 0.5, (2, 1)).opts(line_dash='dashed', line_width=2, alpha=0.5))

Importante

Para poder recuperar perfectamente el espectro continuo (azul) a partir del espectro discreto (rojo) necesitamos que todas sus frecuencias distintas de cero estén dentro del rango \([-F_s/2, F_s/2]\) (negro punteado)

Esta es la base del siguiente teorema fundamental

4.2. Teorema del muestreo

El siguiente resultado se conoce como el Teorema del muestreo

Sea una señal continua \(s(t)\) muestreada a \(F_s\) [Hz] produciendo una señal digital \(s[n] = s(t = n/F_s)\)

La señal continua puede ser recuperada sin pérdidas a partir de las muestras digitales siempre y cuando

\[ f_{\text{max}} < \frac{F_s}{2}, \]

donde \(f_{\text{max}}\) es la componente de frecuencia más alta presente en \(s(t)\)

La recuperación de la señal continua está dada por

\[ s(t) = \sum_{n=-\infty}^{\infty} s[n] \text{sinc}(\pi F_s (t - n /F_s) ) \]

La frecuencia

\[ F_{Nyq} = \frac{F_s}{2} \]

se conoce como Frecuencia de Nyquist, en honor a Harry Nyquist

4.2.1. Ejemplo: Espectro de una señal Gaussiana discreta, parte 2

Asumamos que la frecuencia de muestreo se mantiene en \(F_s=2\) [Hz]. La frecuencia máxima de la gaussiana es la “última” frecuencia donde el espectro es distinto de cero

A continuación se muestra la misma figura anterior pero marcando con un círculo negro la posición aproximada de la frecuencia máxima.

espectro_continuo = hv.HoloMap(kdims=['sigma'], sort=False)
espectro_discreto = hv.HoloMap(kdims=['sigma'], sort=False)
fmax_plot = hv.HoloMap(kdims=['sigma'], sort=False) 

for sigma in [1, 1/np.sqrt(2*np.pi), 0.25]:
    espectro_continuo[sigma] = hv.Curve((x, gauss_freq(x, sigma)), 
                                        'Frecuencia [Hz]', 'Espectro', label='Continuo').opts(line_width=4)
    espectro_discreto[sigma] = hv.Curve((x, periodic_gauss_freq(x, sigma)), 
                                        'Frecuencia [Hz]', label='Discreto')
    fmax = x[gauss_freq(x, sigma)>1e-2][-1]
    fmax_plot[sigma] = hv.Scatter(([-fmax, fmax], [0, 0])).opts(size=10, color='k')
    
(espectro_continuo * espectro_discreto * fmax_plot * hv.Box(0, 0.5, (2, 1)).opts(line_dash='dashed', line_width=2, alpha=0.5))

Cuando la frecuencia máxima sale del cuadrado negro se produce un “solapamiento” de las gaussianas en el espectro discreto

Nota

Debido al solapamiento en el espectro discreto (rojo) se vuelve imposible recuperar el espectro original (azul) sin alteraciones.

El solapamiento espectral se llama aliasing

4.3. Aliasing

Como vimos anteriormente el espectro de una señal muestreada es periódico en \(F_s\). Esto significa que si originalmente la señal tenía componentes con frecuencias mayores a \(\frac{F_s}{2}\) se produce un “traslape” o “solapamiento” espectral

Este fenómeno se llama aliasing y los componentes traslapados se denominan aliases

Importante

Cuando existe aliasing veremos que no es posible reconstruir la señal original sin ambiguedad

4.3.1. Aliasing en el espectro de una sinusoide

Consideremos por ejemplo la siguiente señal sinusoidal

\[ s(t) = \cos(2\pi f_0 t) \]

La transformada de Fourier de coseno es un impulso en \(f_0\) y otro en \(-f_0\)

\[ S(f) = \frac{1}{2} \left(\delta(f-f_0) + \delta(f+f_0) \right) \]

Puedes ver la demostración de esta transformada aquí

Si muestreamos la señal a una frecuencia \(F_s\) que sea menor a \(2 f_0\) entonces habrá traslape en el espectro discreto. Por ejemplo consideremos \(f_0 = 1.23\) [Hz] y \(F_s = 2\) [Hz]. Digamos además que observamos la señal por \(100\) [s]

¿Cómo se ve el espectro de amplitud?

import scipy.fft as sfft

f0 = 1.23 # Frecuencia de la sinusoide
T = 100 # Largo temporal
Fs = 2 # Frecuencia de muestreo

time = np.arange(0, T, step=1/Fs)
signal = np.cos(2.0*np.pi*f0*time)
SA = sfft.fftshift(np.absolute(sfft.fft(signal)))
freq = sfft.fftshift(sfft.fftfreq(n=len(time), d=1/Fs))
hv.Curve((freq, SA), 'Frecuencia [Hz]', 'Espectro de amplitud')

Nota

La señal era de \(1.23\) [Hz], sin embargo en el espectro aparecen peaks en \(\pm 0.77\) [Hz]

Recordemos que el espectro de la señal de la señal continua se repite cada \(F_s\) [Hz] debido al muestreo. La siguiente figura muestra como se ve el espectro discreto de esta señal

../../_images/aliasing2.png

Lo que observamos al calcular la FFT está en el rango

\[ \left[-\frac{F_s}{2}, \frac{F_s}{2} \right] \]

que en este caso es \([-1,1]\) Hz, es decir que no estamos viendo el espectro original (verde) sino un alias (rojo)

4.3.2. Frecuencia de los aliases

Hasta ahora no hemos explicado porqué el alias en este caso tiene una frecuencia de \(\pm 0.77\) Hz

Del diagrama anterior podemos ver que para una sinusoide con frecuencia fundamental \(f_0> \frac{F_s}{2}\), su espectro discreto tendrá infinitos aliases con las siguientes frecuencias

  • \(f_m = f_0 + m F_s\) [Hz] donde \(m\) es un número natural

  • \(f_m = m F_s - f_0\) [Hz] donde \(m\) es un número natural

Según esta fórmula el “primer” alias positivo tiene una frecuencia de

\[ f_m = F_s - f_0 = 2 - 1.23 = 0.773 [Hz], \]

que es exactamente la frecuencia que vimos cuando calculamos la FFT de la señal

4.3.3. Aliasing y reconstrucción de señales

La siguiente figura muestra simula el alias que ocurre cuando muestreamos la señal sinusoidal con frecuencia \(f_0=1.23\) [Hz] usando una frecuencia de muestreo de \(F_s=2\) [Hz]

f0 = 1.23 
T = 5 
Fs = 2 

t_c = np.arange(0, T, step=1e-4)
t_d = np.arange(0, T, step=1/Fs)
s_c = np.cos(2.0*np.pi*f0*t_c)
s_d = np.cos(2.0*np.pi*f0*t_d)
s_a = np.cos(2.0*np.pi*(Fs-f0)*t_c)
p_continuo = hv.Curve((t_c, s_c), 'Tiempo [s]', 'Señal', label='Original').opts(color='g')
p_discreto = hv.Scatter((t_d, s_d), 'Tiempo [s]', label='Muestreada').opts(size=10, color='k')
p_alias = hv.Curve((t_c, s_a), 'Tiempo [s]', label='Alias').opts(line_dash='dashed', color='r')
(p_continuo * p_alias * p_discreto)

Si sólo tenemos la señal muestreada (puntos negros) no es posible distinguir si la sinusoide real es la original (verde) o el alias (roja punteada)

Importante

Debido al aliasing ya no podemos reconstruir la señal sin ambigüedad

4.3.4. ¿Cómo eliminamos el aliasing?

Para que no haya aliasing necesitamos que todas las componentes frecuenciales con amplitud distinta de cero de la señal cumplan con \( |f| < \frac{Fs}{2}\)

Luego, para eliminar el aliasing podemos

  • Modificar \(F_s\): Podemos aumentar \(F_s\) tal que sea dos veces mayor que la frecuencia máxima de interés. Para esto necesitamos saber a priori cuál es la frecuencia máxima

  • Filtrar: Podemos eliminar las frecuencias mayores a \(\frac{F_s}{2}\) antes de muestrear

En la próxima unidad veremos en detalle como filtrar señales

4.4. Apéndice: Principio de incertidumbre para señales

Anteriormente vimos que el ancho en frecuencia y en tiempo de la gaussina son inversamentes proporcionales. Pero en realidad esto es algo que se cumple de forma general para todas las señales.

Recordemos primero el principio de incertudimbre “original”, es decir el de mecánica cuántica

El principio de incertidumbre de Heisenberg nos dice que la precisión (certeza) con que medimos la posición de una particula es inversamente proporcional a la precisión con que medimos su momentum lineal:

\[ \Delta x \Delta p \geq \frac{h}{4\pi}, \]

donde \(h\) es la constante de Planck

Nota

En señales existe un principio análogo: No podemos especificar con infinita precisión la localización temporal y frecuencial de una señal al mismo tiempo.

Denis Gabor (1946) fue el primero en darse cuenta de que el principio de incertidumbre aplica para señales.

Su teorema dice que para una señal con energía finita

\[ E = \int |s(t)|^2 dt \]

con valor medio temporal

\[ \langle t \rangle = \frac{1}{E} \int t |s(t)|^2 dt, \]

y varianza (ancho) temporal

\[ (\Delta t)^2 = \frac{1}{E} \int (t - \langle t \rangle)^2 |s(t)|^2 dt, \]

cuya transformada de Fourier \(\mathbb{FT}[s(t)] = S(\omega)\) tiene un valor medio en frecuencia

\[ \langle \omega \rangle = \frac{1}{E} \int (\omega - \langle \omega \rangle) |S (\omega)|^2 d \omega \]

y varianza frecuencial

\[ (\Delta \omega)^2 = \frac{1}{E} \int (\omega - \langle \omega \rangle)^2 |S(\omega)|^2 d\omega \]

se cumple que

\[ \Delta t \Delta \omega \geq \frac{1}{2}, \]

es decir \(\Delta t\) y \(\Delta \omega\) no pueden ser arbitrariamente pequeños

Importante

El ancho temporal y el ancho frecuencial están inversamente correlacionados sin importar la señal en particular